
__device__ prec computeDistances(float3 ri, float3 rj) { 

	prec r, cp, cx, cy, cz; 

	cx = ri.x - rj.x;
	cy = ri.y - rj.y;
	cz = ri.z - rj.z; 

	//apply periodic boundary conditions
	cx = cx - box.x*round(cx/box.x);
	cy = cy - box.y*round(cy/box.y);
	cz = cz - box.z*round(cz/box.z);

	cp = cx*cx + cy*cy + cz*cz;

	r = sqrt(cp);

	return r; 
}


__host__ prec radius(int i, int j, int k, int l){

	prec r, cp, cx, cy, cz;

	//get r
	cx = (hAtom[i*3+k].x-hAtom[j*3+l].x);
	cy = (hAtom[i*3+k].y-hAtom[j*3+l].y);
	cz = (hAtom[i*3+k].z-hAtom[j*3+l].z);

	//apply periodic boundary conditions
	cx = cx -box.x*round(cx/box.x);
	cy = cy -box.y*round(cy/box.y);
	cz = cz -box.z*round(cz/box.z);

	cp=cx*cx+cy*cy+cz*cz;
	r = sqrt(cp);

	return r;

}

__host__ prec cpuEnergy(double r, int i, int j, int k, int l){

	prec sig, eps;
	prec dP = 0.0;
	prec a = 2.4829;//cutoff term
	prec B = 2.7956455;

	if(k !=2 && l == 2 || k == 2 && l != 2) {
		sig = sigoc;
		eps = epsoc;
	}
	if(k != 2 && l != 2) {
		sig = sigoo;
		eps = epsoo;
	}
	if(k == 2 && l == 2) {
		sig = sigcc;
		eps = epscc;
	}

	prec rcut = a*sig;

	if(r < rcut) {
		prec f1 = sig/r;
		prec f2 = f1*f1;
		prec f6 = f2*f2*f2;
		prec xp = exp(sig/(r-a*sig));
		prec f12 = B*f6*f6;
		dP = eps*((B*f6*f6)-f6)*xp;
	}
	else dP = 0;

	return dP;
}

__device__ prec gpuEnergy(prec r) { 

	prec dP = 0.0;
	prec a = 2.4829;//cutoff term
	prec B = 2.7956455;

	prec RCUT = a*sig;

	if(r < RCUT) {

		prec f1 = sig/r;
		prec f2 = f1*f1;
		prec f6 = f2*f2*f2;
		prec xp = exp(sig/(r-a*sig));
		prec f12 = B*f6*f6;
		dP = eps*((B*f6*f6)-f6)*xp;
	}
	else dP = 0;

	return dP;
}

__device__ prec phi0(float3* hAtom) {

	prec sig, eps; 
	prec P1 = 0.0;
	prec dP,r;	

	int idx = blockIdx.x * blockDim.x + threadIdx.x; 

	//loop to choose molecule and atom ID 
	for (int j=0; j < blockDim.x; j++) { 
		for(k=0; k<=2; k++) {
			float3 ri = hAtom[idx*numSite+k];
			for(l=0; l<=2; l++) {
				if (k!=2 && l==2 || k==2 && l!=2) { sig = sigoc; eps = epsoc; }
				if (k!=2 && l!=2) 		    { sig = sigoo; eps = epsoo; }
				if (k==2 && l==2) 		    { sig = sigcc; eps = epscc; }

				float3 rj = sharedPos[threadIdx.x*numSite+l];	
				r = computeDistances(ri, rj);
				dP = gpuEnergy(r, sig, eps);
				P1 = P1 + dP;
			}
		}	
	}

	return P1;
}

__global__ void CO2Energy(void *rAtom, int numMolecules, int ref, float3 box) { 
	
	
	int idx = blockIdx.x * blockDim.x + threadIdx.x; 

	for (int i=0, tile=0; i < numMolecules; i+=p, tile++) { 
		
		sharedPos[threadIdx.x*numSite+0] = gAtom[idx*numSite+0];   // Carbon
		sharedPos[threadIdx.x*numSite+1] = gAtom[idx*numSite+1]; // Oxygen
		sharedPos[threadIdx.x*numSite+2] = gAtom[idx*numSite+2]; // Oxygen
		
		__syncthreads(); 

		if (ref==0) {  
			twobodyEnergy += phi0(); 
		} else {
			twobodyEnergy = phi1(moveSelect);
		}
		__syncthreads(); 

	}
}
		
__host__ prec phi1(int moveSelect) {
	
	prec P1 = 0.0;
	prec dP,r;

	int i = moveSelect;

	for (int j = 0; j<=(nMolecules -1); j++) {
		if (j != i) {
			for(k = 0; k<=2; k++) {
				for(l = 0; l<=2; l++) {
					r = radius(i, j, k, l);
					dP =  cpuEnergy(r,i, j,  k, l);
					P1 = P1 + dP;
				}	
			}
		}
	}
	
	return P1;
}

prec calculatePotential(float3 *hAtom) {

	// Calculate change in energy on the CPU
	if(ref == 1) {
		pot = phi1(moveSelect);
	} else {	// Compute total energy on the GPU 
				// HOST --> DEVICE MEMORY TRANSFER 

		cudaMemcpy(dAtom, hAtom, nMolecules*3*sizeof(float3), cudaMemcpyHostToDevice);

		int blockDim = (nWaters/GPUTHREADCOUNT)+1;
		printf("Number of blocks = %d\n", blockDim);

		cyclic_EnergyCalculation<<<blockDim, GPUTHREADCOUNT>>> (dAtom, d2body, nMolecules, box);
		cudaThreadSynchronize();

		// DEVICE --> HOST MEMORY TRANSFER

		cudaMemcpy(h2body, d2body, nWaters*sizeof(precd), cudaMemcpyDeviceToHost);

		precd gpuEnergy_1 = 0.0;

		for (int i=0; i<nWaters; i++)  {
			gpuEnergy_1 += h2body[i];
		}
		
		pot = gpuEnergy_1; 
	}
	printf("moveSelect = %d, ref = %d, Energy = %g\n", moveSelect, ref, pot); 

	return pot*ENERGY_CF;	
}


